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ABSTRACT 

We analyze observations of the Seyfert galaxy NGC 4151 covering 90 years in the 
optical band and 27 years in the 2-10 keV X-ray band. We compute the Normalized 
Power Spectrum Density (NPSD), the Structure Function (SF) and the Autocorre- 
lation Function (ACF) for these data. The results show that the optical and X-ray 
variability properties are significantly different. X-ray variations are predominantly in 
the timescale range of 5 - 1000 days. The optical variations have also a short timescale 
component which may be related to X-ray variability but the dominant effect is the 
long timescale variability, with timescales longer than ~ 10 years. We compare our re- 
sults with observations of NGC 5548 and Cyg X-1. We conclude that the long timescale 
variability may be caused by radiation pressure instability in the accretion disk, al- 
though the observed timescale in NGC 4151 is by a factor of few longer than expected. 
X-ray variability of this source is very similar to what is observed in Cyg X-1 but scaled 
with the mass of the black hole, which suggests that the radiation pressure instability 
does not affect considerably the X-ray production. 

Key words: galaxies: active - galaxies:individual:NGC 4151 - instabilities - X- 
raysrgalaxies. 



■ 1 INTRODUCTION 



NGC 4151 {z = 0.00332; D = 13.2 Mpc assuming Hubble 
constant 75 km s~^ Mpc~^) is one of the classical Seyfert 
1 galaxies. It is among the most frequently observed AGN 
and therefore it is an excellent object for variability stud- 
ies, although it is both typical and atypical representative 
of its class, as nicely and extensively summarized by Ul- 
rich (2000). We therefore chose it as a second source (after 
NGC 5548; Czerny, Schwarzenberg-Czerny, & Loska 1999) 
for which the nature of the variability can be studied by 
computing the normalized power spectrum density (here- 
after NPSD) in both the optical and the X-ray bands. 

The host galaxy is a typical spiral galaxy SAB(rs)ab 
of 11.3 mag brightness in B band (Perez et al. 1998). The 
extinction due to our Galaxy is low {Ay < 0.01 in B band) 
because NGC 4151 is located near the North Galactic Pole 
(6// — 75.1°). However, the intrinsic extinction in NGC 4151 
is considerable: in the optical region, Av ~ 1™.0 in BLR and 
Ai^ = 0^.7 in NLR (Ward et al. 1987). Direct estimates in 
the UV bands give much less extinction (Eb-v ~ 0.06) for 
wavelength A2200 A corresponding to Av = 0'".18 (Wu 



& Weedman 1978). The hydrogen column density inferred 
from the equivalent widths of ultraviolet absorption CHI 
lines is about 1.7 x 10^'"* cm~^ (Kriss et al. 1992), smaller 
than the inferred X-ray column density. The X-ray obser- 
vations clearly indicate that the nuclear X-ray continuum 
is significantly obscured, by an equivalent hydrogen column 
Nh = {1- 10) X 10^^ cm"^) (Yaqoob et al. 1993), which 
would correspond to Av ~ 5™ — 50™ according to the stan- 
dard relation between Av and Nh from Reina & Tarenghi 
(1973). 

HST observations suggested that the line of sight to- 
wards the nucleus (at least in the epoch when the observa- 
tion was performed, on 1991 June 18/19) lies outside the 
ionization cone (Evans et al. 1993). However, broad emis- 
sion lines - which can vary in strength - are usually eas- 
ily detected, although with superimposed absorption fea- 
tures (e.g. Brandt et al. 2001); only during one observation, 
in 1984, the broad component of H/j was absent (Lyutyj, 
Oknyanskij & Chuvaev 1984; Penston & Perez 1984; see 
also Pronik, Sergeev, & Sergeeva 2001) and the spectrum 
more closely resembled that of a narrow emission line or a 
Seyfert 2 type galaxy. Therefore it seems that most of time 
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the nucleus is viewed directly and the scattered component 
characteristic for Scyfcrt 2 galaxies does not dominate the 
spectrum, apart from the soft X-ray band (see Ogle et al. 
2000). 

Regardless of the fact that X-ray extinction is high, 
the intrinsic X-ray spectrum seems to be fairly typical for 
a Seyfert galaxy (Zdziarski, Johnson, & Magdziarz 1996). 
The brightness of the nucleus varies dramatically on all 
timescales in all energy bands. Spectrum and variability in 
the 50 - 150 keV band was studied by Johnson et al. (1997); 
the source shows a factor of 2 variability in days and years. 
The X ray (2 - 10 keV) variability of this source was studied 
in more detail using Fourier analysis, starting with the paper 
by Fiore et al. (1989). It was found (Papadakis & McHardy 
1995) that X-ray variability power spectrum of NGC 4151 
covering the range 10"^ - 10"* Hz can be described by a 
power law, which is significantly flatter than that measured 
on short time scales (10"^ -10"^ Hz) derived from EXOSAT 
observations. The power law index measured by Hayashida 
et al. (1998) on the basis of the Ginga data in the frequency 
region of lO"** — 10~^ Hz has an index a = 2.1. 

The long-term optical variability of the continuum in 
the UBV band was studied by Lyuty & Doroshenko (1999). 
They showed that a long photometric minimum in years 
1984 - 1989 separates two active phases: 1968 - 1988 (Cycle 
A) and 1989 - 2000 (Cycle B). The decrease of brightness 
continued in 2000 and at the end of 2000, the brightness of 
the variable component reached almost the same level as in 
1987 - 1988 (Doroshenko et al. 2001). These phases are also 
seen in the JHKL-bands (Lyutyi, Taranova, & Shenavrin 
1998). 

The power spectrum density (PSD) of the NGC 4151 
optical flux variations in 1968 - 1987 was studied by Tere- 
bizh, Terebizh, & Biryukov (1989). It was found that the 
PSD in the frequency range 10~* — 2 X 10~^ day~^ (time in- 
terval from 50 days to 30 years) corresponds to flicker noise 
with slope « 1 in logarithmic scales and the light curve may 
be interpreted as a result of superposition of flares randomly 
distributed in time. The behavior of the structure function 
(hereafter SF), focusing on the intra- night variability and 
conducted between 1989 - 1996 was studied by Merkulova, 
Metik & Pronik (2001). Collier & Peterson (2001), using the 
same technique, determined the characteristic timescale to 
be 13t" days. 

Long timescale trends in optical variability of NGC 4151 
were studied by Lyutyj & Oknyanskij (1987) and Fan & 
Su (1999), on the basis of data covering the period from 
1910. Lyutyj & Oknyanskij (1987) noted the quasi-periodic 
changes with typical time of 4 and 14 years. Also Fan & 
Su (1999) claimed the presence of the periodicity, with the 
period 14.0 ± 0.8 yr, during which both the brightness and 
the spectral type changed. It is not surprising that both 
groups found the similar quasi-periodical processes in optical 
light curve because the main data set used by both teams 
was the same. Some earlier papers also claimed the presence 
of the periodic variability albeit with other values of the 
period (e.g. Pacholczyk (1971) found a period of 5.1 years), 
while other papers argued against any strict periodicity (e.g. 
Pacholczyk et al. 1983; Terebizh et al. 1989). 

The short term multiwavelength variability was best 
studied by AGN Watch team and was summarized by Edel- 
son et al. (1996). They concluded that significant and corre- 



lated variability was observed in the X-ray, UV and optical 
bands, with phase differences consistent with zero lag and 
normalized variability amplitude decrcEising with increasing 
wavelength. They calculated the power spectrum in the UV 
on timescales of « 0.2 - 5 days and showed that the power 
spectrum is falling down rapidly at short timescales and 
the bulk of the variability power is on timescales of days 
or longer. 

In the present paper we reanalyze the optical and X-ray 

data available in the literature which cover respectively 90 
and 25 years. Such a broad dynamical range allows us to 
draw conclusions about the character of the accretion fiow 
and the possible nature of the instabilities responsible for 
the observed variability. 



2 OBSERVATIONAL DATA 
2.1 Optical band 

The observational data in optical region used in the present 
paper come from different sources. Long data sets were col- 
lected by the team at the Crimean Station of the Stern- 
berg Astronomical Institute (1968 - 2000), the Special As- 
trophysical Observatory in the Caucasus (1997 - 2000), and 
the Maidanak Observatory of the Ulugbek Astronomical 
Institute in Uzbekistan (1990 - 2000). Below, we will re- 
fer to these data as the Crimean data. The data collected 
in the three bands (U, B, and V) cover the period from 
1968 till 2000. The data points are not distributed uni- 
formly since the source was at some seasons unobservable 
at the sites involved. Nonetheless, the sets consist of about 
1150 points in each color, thus providing an exceptionally 
long and well sampled light curve for any AGN. The light 
curve and the discussion of the variability amplitude and 
color changes were presented by Lyuty & Doroshenko (1999) 
and Doroshenko et al. (2001). Additionally, shorter data set 
(about 300 points during 1989 - 1999) in the infrared band 
R is also available (Doroshenko et al. 2001). 

Shorter, but more densely spaced data set was taken 
from the AGN Watch team (Kaspi et al. 1996). Those data 
cover only 98 days but they add 62 points to the long set. 
Continuum flux in the spectral band 4560 - 4640 A was fitted 
to the flux in the B-band of Crimean data set using regres- 
sion relation based on 17 common dates of observations. 

The majority of the oldest data points was derived by 
Pacholczyk et al. (1983) from the old Harvard and Stew- 
art Observatories photographic plates. These data cover the 
period from 1910 until 1968. Besides those observations we 
used data from Fitch, Pacholczyk, & Weymann (1967), Can- 
non, Penston, & Brett (1971), Pension et al. (1971), and 
Oknyanskij (1978; 1983). Comparison of the photographic 
magnitudes between various data sets and with Crimean 
photoelectric data in B-band for overlapping dates of ob- 
servations allowed us to reduce the photographic data into 
photoelectric B-bands magnitudes. 

As a result, in the B-band we obtained a quite long data 
set covering years 1910 - 2000. We must note that the similar 
work of gathering all photographic observations was made by 
Lyuty and Oknyanskij (1987), but in comparison with this 
work we added more photoelectric observations made after 
1984. A subsequent compilation of the 1910 - 1991 data was 
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prepared by Longo et al. (1996) but they used even older 
data from heterogeneous sources and interpolated the data 
in order to obtain equally spaced data points which lead to 
large errors so finally they had to ignore the first 30 years 
of observations. 

The full combined light curve of the NGC 4151 nucleus 
in the B band is shown in Fig. Open circles mark the his- 
torical photographic plate data while the Crimean data are 
shown with the continuous line. A small expanded fraction 
of the light curve from JD2449150 till JD2449550 is inserted 
in that Figure. It includes the three month period when the 
dense AGN Watch data were collected. Crimean data are 
shown with continuous line and the AGN Watch points are 
marked with crosses. The figure demonstrates that the vari- 
ability in 1910 - 1978 based on photographic data is similar 
in its character to the variability in 1968 - 1989 covered 
by photoelectric observations. Also the AGN Watch data 
agrees well with the corresponding fragment of the Crimean 
observations: the AGN Watch data points are located on 
the quasi-linear part of the light curve showing local rise in 
brightness. Therefore, we consider the old data as generally 
reliable. 

The overall character of the variability is not constant 
over the entire data set - the behavior of the optical light 
curve in 1990 - 2000 appears different from that in 1910 - 
1985. The evolution up to 1985 can be represented by a sin- 
gle parabolic-type curve with maximum fiux of about 75 - 80 
mjy in B band and with superimposed on it more rapid vari- 
ations. In 1935 there was (if the photographic magnitudes 
were estimated correctly) the strongest flare which lasted 
only 10 - 15 days. A less intense flare happened in 1946 but 
it is based on a single point. The evolution in 1990 - 2000 
appears as a single giant outburst but this change in the 
character of variability occurred in the middle of the epoch 
1968 - 2000 covered by the photoelectric data. 
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Figure 1. The light curve of NGC 4151 in the B band from the 
historical plate data (open circles connected with a continuous 
line) and Crimean observations (solid line). Inserted expanded 
fragment shows the AGN Watch data (crosses) together with 
Crimean data (continuous line). The errors are on the order of 
0.2 mag (~ 20 %) for the historical plate data and ~ 1.5 % for 
the Crimean observations. 
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2.2 X-ray data 

The data in the 2-10 keV X-ray band also came from 
different papers and databases which were collected using 
different satellites. Fluxes were corrected for absorption. The 
observed time coverage is from October 1974. 

The data points during 1975 - 1992 were taken from 
Papadakis & McHardy (1995). Of their 133 data points, a 
majority came from the Ariel V Sky Survey, and others were 
taken from literature including observations by EXOSAT, 
GINGA, OSO-8, HEAO-1, TENMA and Ariel VI satellites. 
Some additional points (in period May 1987 - May 1995) 
were taken from Yaqoob & Warwick (1991) and Yaqoob et 
al. (1993) and were also taken from the Tartarus Database 
of ASCA observations ^. We corrected their fluxes for ab- 
sorption using the HEASARC's onhne W3PIMMS Version 
3.1 flux converter ^ as well as papers of Weaver et al. (1994) 
and Edelson et al. (1996). Data from the first three years of 
RXTE observations during October 1998 - November 1999 
were taken from Markowitz & Edelson (2001). The observa- 
tional data given in that paper in counts/ s were transformed 
into unabsorbed fluxes via the mean luminosity of 300-day 



^ http://tartarus.gsfc.nasa.gov/ 
^ http://hcasarc.gsfc.nasa.gov/ 
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Figure 2. The light curve of NGC 4151 (vFy) in U band (solid 
line) and in X-ray band (2 - 10 keV absorption-corrected flux, 
open circles). The errors of X-ray data are on the order of a 
factor 2 in the oldest data points and of ~ 10 % for more recent 
measurements. The errors in Crimean data are ^ 2.2 %. 



window given in the paper. We supplemented those gathered 
data with three long sequences. Two EXOSAT light curves 
were obtained from Ian M. McHardy for the paper by Cz- 
erny & Lehto (1997). These data cover 10 and 11 July 1983 
and 15 May 1985, with rate sampling every 100 s. Third 
sequence, from the ASCA satellite, covers the period from 
12 till 23 May 2000, and is nearly continuous (except for 
Earth occultations and the passages through the South At- 
lantic Anomaly). For that ASCA data set, we use the count 
rate sampling of 32 s. The measurements in counts/s were 
again converted to fluxes using the mean luminosity for an 
appropriate data sequence. 

All daily-averaged X-ray points used by us for further 
analysis are shown in Fig. |21 with open circles. 
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3 METHODS 



3.2 Structure function 



3.1 Power spectrum 

We attempt the restoration of the underlying power spec- 
trum |_F|^ of the active nucleus from its observed light 
curves, i.e. from the product of the true light curve and the 
sampling (window) function. Direct analysis of observations 
yields the observed power spectrum |Gp, i.e. the underly- 
ing spectrum affected by the window function These 
functions obey an exact relation in the Fourier transform 
space: 

G = F*W (1) 

where * denotes the convolution and F, G and W are Fourier 
transforms of the underlying and observed light curves and 
of the sampling function. For long and nearly uniformly sam- 
pled observations \W\^ approaches the delta function and G 
approximates F. For non-uniform sampling, the shape of 
the window function becomes complex, with broad wings 
and many local maxima. In such case it is difficult to obtain 
a solution of Eq. Q for F as it involves deconvolution of 
the sampling transform W . 

We address the issue of uneven sampling by not cal- 
culating the power spectrum of the entire light curve di- 
rectly. Instead, we analyze the properties in various time 
scale ranges separately. 

In the optical band, we first rebin the data to 1 yr bins 
and fill the few existing gaps by linear interpolation, thus 
obtaining 90 data points. Next we rebin the original curve 
into 5 day bins and find an almost continuous sequence con- 
taining 2403 such points. Data gaps are again filled through 
interpolation. We calculate the power spectra of the two 
sequences separately and combine them by averaging the 
logarithmic values. 

In the X-ray band, we also first rebin the data to 1 yr 
bins (25 points). Next we rebin the original light curve into 
30 day bins and find and almost continuous fragment (26 
points). Next we rebin the original data into 5 day bins ob- 
taining 69 points in an almost continuous set. Finally, the 
two light curves from EXOSAT and one from the ASCA 
monitoring were rebinned to 4096 s bins. The power spec- 
tra of the resulting six light curves were computed indepen- 
dently and combined for further analysis. 

For uniformity, all optical data points are converted 
from magnitudes to flux units, when necessary, and the av- 
erage value is subtracted. The NPSD of each light curve is 
computed according to the description by Hayashida et al. 
(1998): 



NPSD{f)^=^a{f)a{f)xT, 
x{t) 

with 

1 " 

«(/■) = - x{ti) exp{2nifti), 



(2) 



(3) 



where x{ti) is the value of the flux at the time ti, T is the 
length of the data set and x{t) is the mean value of the fiux. 
Therefore, the integrated NPSD gives half of the fractional 
variance. 

The errors of the power spectra obtained as above are 
either estimated directly from the data or estimated using 
the Monte Carlo simulations (see Appendix A). 



For the analysis of our time series we also used the structure 
function technique (SF). The application of SF allows us 
to study both stationary and nonstationary processes. The 
SF was introduced by Kolmogorov in 1941 for analysis of 
the statistical problems connected with turbulence theory 
(Kolmogorov 1941a; b). In astronomical practice, the wide 
application of the SF to time-series analysis begun in the 
middle eighties (e.g. Simonetti, Cordes, & Heeschen 1985; 
Paltani, Courvoisier, & Walter 1998; Kataoka et al. 2001). 

As a rule, most work in astronomy uses only the first- 
order structure function (SFi). By definition the structure 
function of the first-order is a mean square of difference 
[x{t) ~ x{t + r)]: 



SFi{r) = M[{x{t) - x(t + r) 



(4) 



where x{t) is random process and r is time shift. Since we do 
not use higher order structure functions, we simply denote 
SFi as SF. 

In the case of a stationary random process the SF is 
directly related to the autocorrelation function: 



SF{t) = 2D[x{t)] X [1 - ACF(r) 



(5) 



where D[x(t)] is variance of the process and ACF is the 
autocorrelation function. 

The slope of the SF changes with time interval t. If the 
measurement errors are neglected, on the shortest time scale 
the variability can be well approximated by a linear trend, 
and then SF cc r^. For long time scales, the slope of the SF 
becomes fiatter and in the limit when r — » oo the structure 
function saturates: SF — > 2D[x{t)]. The addition of mea- 
surement errors to the random process increases SF by the 
value 2Derr where Derr is the variance of the measurement 
errors. Therefore, SF 2Derr when r — > 0. 

Many processes are stochastic in nature. For some pro- 
cesses such as fractional Brownian motion, the SF shows a 
power law dependence on r. Some processes can be repre- 
sented as a superposition of a large number of random im- 
pulses of deterministic form (shot noise). Sub-group of those 
processes which are characterized by the PSD{f) oc l/f 
are called flicker noise. If 7 takes values from 1 to 3, then in 
this case the structure function SF{t) oc where b takes 
values from to 2, following the relation 



7 = 6-1-1. 



(6) 



The structure function analysis was done using a soft- 
ware package by S.G. Sergeev from the Crimean Astrophys- 
ical Observatory. The whole time interval was divided into 
equal bins in logarithmic scale and for each bin we found 
such pairs of observations with tj > ti that their time differ- 
ence Tfc — tj — ti fltted into the given bin. Next we calculated 
the value 



SF{r,) = 



[x{t,) - x{t,)]^ 
rik 



(7) 



where nt is the number of pairs in k-th bin. We estimate 
the errors of the SF as described in the Appendix B. 

To allow a meaningful comparison of the SF in various 
bands we introduce also the normalized structure function 
(NSF) as follows: 
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Table 1. The basic parameters of the variabiUty. 



Color band Flux (mjy) rms (mjy) Fvar R max error 

Whole data 

U,B,V (1968-2000), R (1989-2000), X-ray (1975-2001), Btot (1910-2000) 



X-ray 


[2-10] keV 


2.55e-2* 


1.07e-2* 


0.41 


37.4 


10.0% 


U 


3600 A 


48.44 


24.25 


0.50 


8.7 


2.5% 


B 


4400 A 


63.56 


21.09 


0.33 


3.9 


1.5% 


Btot 


4400 A 


63.23 


19.17 


0.32 


5.2 




V 


5500 A 


92.26 


19.45 


0.21 


2.6 


1.0% 


R 


7000 A 


169.84 


30.78 


0.18 


2.5 


2.0% 



Cycle A (25.03.1968 - 03.07.1988) 








U 3600 A 31.88 


10.66 


0.33 


2.5 % 


B 4400 A 48.26 


9.10 


0.18 


1.5 % 


V 5500 A 77.82 


8.37 


0.11 


1.0 % 


Cycle B (11.02.1989 - 13.03.2000) 








U 3600 A 63.51 


23.32 


0.37 


1.3 % 


B 4400 A 75.35 


20.09 


0.27 


0.74% 


V 5500 A 103.43 


18.18 


0.18 


0.83 % 


Fvar - rms to mean ratio, Rmax 


- maximum 


to minimum 


flux ratio 


* Flux at 2 keV computed assuming 


a = 1 energ 


;y index in 2 - 


■ 10 keV band 



In this case all NSF should approach the universal value 2 
on long timescales. 

3.3 Autocorrelation function 

In our computation of ACF, we again use of program pack- 
age by S.G. Sergeev; the method used relies on the inter- 
polation. Specifically, for any point from the real data, a 
corresponding shifted point is found by linear interpola- 
tion or extrapolation, and extrapolation is made adapting 

x(^t > tmax) ~ xi^tmax) and x[t < train) ~ x(^tmin) (for a 

more detailed description, see Sergeev et al. 1999). 



in X-ray band while the systematic long timescale dim- 
ming/brightening seems to be stronger in U band than in 
X-rays. The energy content is significantly greater in the 
variable part of the optical spectrum than in the X-ray spec- 
trum. Variable part of uFi, in U band is equal to 2.02 x 10"^" 
erg cm~^ s~^ while the variable part of the 2-10 keV fiux 
is only 0.81 x 10^^" erg s~^ cm~^. Of course we do not know 
the bolometric corrections needed for more strict comparison 
of the total X-ray and optical variable energy output but, 
they are not likely to reduce this discrepancy considerably. 
Therefore, we conclude that strong long timescale trends in 
the optical band do not result from reprocessing of the vari- 
able X-ray emission, unless X-ray variability leads to strong 
variations of the extinction, thus amplifying the effect. 



4 RESULTS 

4.1 General variability properties of NGC 4151 

The basic variability properties of the Seyfert galaxy 
NGC 4151 seen in the data assembled by us are summa- 
rized in the upper part of the Table [H We give there the 
mean fiux at the given spectral band, rms fiuctuations, the 
normalized amplitude F^ar (i.e. rms to mean ratio) and the 
ratio i?max of the maximum to the minimum flux. The full 
90 years of the B band data are given separately from the 
data subset determined by photoelectric measurements. 

The normalized variability amplitude is the largest in 
U band and decreases towards the longer wavelengths. This 
effect is well established for AGN (e.g. Ulrich et al. 1997) 
and tells us that the variable component is bluer than the 
non- variable (starlight) component. 

The overall normalized variability amplitude in the X- 
ray band is comparable to the amplitude in U band. How- 
ever, more extreme variability events occasionally happen 



4.2 Stationarity of the data 

Question of whether light curves of AGN covering a time 
span of several decades may be considered realizations of 
a stationary process is interesting per se, but also because 
of its possible influence on investigations of variability on 
shorter time scales (e.g. Leighly 1999, Press & Rybicki 1997). 
We address this issue by estimating of the importance of any 
existing linear trend. 

4-2.1 Optical band 

Even the visual inspection of the light curves of NGC 4151 
clearly shows some long timescale trends although they be- 
come less prominent with an increase of the data length. It 
can be best studied in the B band since at that color the 
available light curve covers the longest time span. 

We analyzed the longest trends by rebinning the data 
in one year bins in order not to be biased by the most re- 
cent time spans with superior data coverage as compared 
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with the earliest data from the beginning of the previous 
century. We studied the trends both using magnitudes (ef- 
fectively logarithmic scale) or using uF,j fluxes (effectively 
linear scale). 

Mean value of the B magnitude in the Crimean obser- 
vations is 11.83. Apart from two major outbursts, the pres- 
ence of the linear trend is also seen in the data. Linear fit 
in magnitude scale gives a shift by -0.35 mag during the 
entire period of observations (32 years). This is comparable 
to the value of the dispersion in this period, calculated for 
unbinned data (0.33 mag). The result is similar if flux {vF^) 
is used: linear trend gives the net brightening by 1.2 x 10"^^° 
erg s~^ cm~^, with the dispersion equal 1.4 x 10"^" erg 
s~^ cm~^. However, if the entire data set of 90 years of 
data is analyzed, the linear trend acts in the same direction 
but is much weaker. Systematic shift is only -0.13 mag, or 
3.6 X 10~^^ erg s~^ cm~^, with the dispersion unchanged. 

This shows that the data are not strictly stationary even 
on the timescale of 90 years. However, in the longest data 
set the linear trend is much weaker than the typical am- 
phtude of the variability (in B band) and the criterion for 
stationarity is roughly satisfied. It is less so in the case of 
shorter sequences, so the power spectra and the structure 
function may depend on the choice of the data set. This 
is particularly well known from the studies of X-ray light 
curves of galactic sources (e.g. Uttley & McHardy 2001, and 
references therein; for quasar sample see Manners, Almaini, 
& Lawrence 2001). 




-4 -2 



log(f) [1/day] 

Figure 3. Normalized Power Spectrum Density of NGC 4151. 
Upper panel: U band (1968 - 2000); middle panel: B band (1910 - 
2000); lower panel: X-ray band (1974 - 2000). Marked errors are 
direct observational errors, as described in Appendix A. 



4.3 Power spectrum 



4.2.2 X-ray hand 

X-ray data cover only the period from 1974 till 2000, and the 
older data (till ~ 1986) seem to show larger scatter so the 
presence or the absence of the long trends is less apparent 
(see Fig. 13. 

We again rebinned the data in one year bins in order to 
achieve a uniform distribution on long time scales and looked 
for a linear trend in the resulting light curve. The mean 2-10 
keV fiux was 1.98x 10"^" erg s~^ cm~^ while the linear trend 
gave the systematic brightening of the source by 4.12 x 10~^^ 
erg s~^ cm~^, much lower than the dispersion in the data 
during this period (8.1 x 10~^^ erg s~^ cm~^). Such absence 
of the linear trend in the data is partially accidental - more 
accurate Ginga/ASCA/RXTE data (1978 - 2000) cover the 
period which started and ended up at a very similar level of 
the activity of the source. Slightly shorter sequences show 
trends up to almost an order of magnitude stronger but still 
smaller than the dispersion. 

This suggests that the X-ray fiux is more likely to be 
produced by a stationary process during the observed period 
than the optical fiux during the same epoch. Note however, 
that the time scales used for sampling of the X-ray fiux range 
from seconds up to decades, i.e. have a greater dynamical 
range than the time scales of optical sampling, which range 
from days up to nearly a century. 

In this context a detection of the change of slope (and 
in particular, of fiattening) of the PSD towards longer time 
scales might be easier in X-rays than it would be in the 
optical band. A few more years of occasional monitoring are 
needed to confirm such a statement. 



4.3.1 Optical band 

We analyze first the Crimean data alone, covering the years 
from 1968 till 2000, in all three color bands. 

The power spectrum in the U band (see Fig. |3 top 
panel) has roughly a power law shape with a steep slope 
(~ 2) but the presence of some structure is clearly visible. A 
flattening is seen at a timescale of ~ 300 d which seems to 
separate two variability components. No flattening is seen 
at the longest timescales but the data set used to determine 
this power spectrum covers only ~ 10 years. In the V band, 
the spectrum (not shown) shows practically no substruc- 
ture. Given all this, we infer that the interesting timescales 
are rather long, and therefore concentrate on the results de- 
rived from the complete set of data in B band, covering the 
entire 90 years. The result is shown in Fig. middle panel. 
In similarity to the U band power spectrum, the B band 
data show a hint of the presence of two components, sepa- 
rated at the timescale on the order of ~ 300 days. The data 
also reveals the flattening of the low frequency component 
(roughly at timescale of ~ 10 years) although the measure- 
ment errors are large. 

In order to determine whether the fiattening of the 
power spectrum is real and not simply caused by the fi- 
nite extension of the data, we performed Monte Carlo sim- 
ulations assuming that the power spectrum of the entire 90 
years of the data in B band is represented by a broken power 
law model 



if log/>/o 
The best fit values of the model parameters are given 



log f < .fo 
log .f > .fo 
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Figure 4. Upper panel: error contours (thick line - 68% confi- 
dence level, thin line - 90% confidence level) for the fit of the 
power law with an index a and a break at log /o to the NPSD of 
NGC 4151 in the B band (1910-2000; dashed line) and X-ray band 
(continuous line) determined from Monte Carlo simulations (see 
Appendix A). Cross marks the location of the best fit. Interme- 
diate and lower panels show the plots of the marginal probability 
for both parameters in the X-ray band (continuous line) and in 
the optical band (dashed line). 



in Table|21 together with 1 a errors for each parameter. Full 
two-dimensional contour errors are shown in Fig. 2] Despite 
the fact that the data used for this cover such an excep- 
tionally long period, the errors are still considerable: the 
NPSD is consistent with a single power law. It is steeper than 
1.5 and the flattening happens at timescales longer than 15 



years. We note here that the outlined procedure for compu- 
tation of the confidence levels is based on likelihood function 
derived from simulations. Conceptually it is different from 
the method employed in Appendix A. The two methods 
converge only in an asymptotic limit of infinite number of 
observations. Here we adopt the former method as the more 
conservative one, yielding much broader confidence limits. 

We reploted the resulting NPSD in the form of Power x 
Frequency diagram (see Fig. |^ in order to better illustrate 
the dominant time scales. The component corresponding to 
the timescale of years appears to dominate, but the Figure 
also shows the presence of a separate, month-timescale com- 
ponent, as well as a possible high frequency roll-over to the 
last component at about log/(days)= —2. It would agree 
with the characteristic timescale determined by Collier & 
Peterson (2001). 

Lyuty & Doroshenko (1999), analyzing the data from 
the period 1968 - 2000, argued that there were two cycles of 
activity in NGC 4151 which difl^ered considerably regarding 
the level of variability: first one from 1968 till 1988 (cycle 
A) and other from 1989 till the present time (cycle B). Since 
the duration of the outburst B - the dominant feature in the 
optical light curve (see Fig. - corresponds to the char- 
acteristic frequency seen in our analysis, we repeated this 
analysis for the data set in the B band without the cycle 
B. The NPSD computed without the cycle B still showed a 
flattening at log/(days)~ —3.6 but above it the NPSD was 
even better represented by a single power law with a slope 
~ 2. 

4.3.2 X-ray band 

X-ray PSD is distinctively different from the optical NPSD 
(see Fig. Most of the power is now on much shorter 
timescales. The spectrum generally shows a flattening on 
time scales of about 100 days, but more structure can be 
identified. The spectrum might be crudely represented as 
a broken power law, with slope 1 at intermediate frequen- 
cies and about 2 at higher frequencies. The normalization 
at high frequencies is similar to that obtained by Hayashida 
et al. (1998). The position of the flattening determined by 
Papadakis & McHardy (1995) is now in the middle of the 
intermediate part. 

The characteristic features of the shape are again ap- 
parent on Power x Frequency plot (see Fig.|SJl. Most of the 
X-ray power is concentrated between ~ 6 days and ~ 150 
days, with a turnover both at higher and lower frequencies. 
There may also be a decrease in the variability level at the 
intermediate timescales of ~ 30 days. There may also be 
some additional power at the lowest frequencies but at the 
lower level than in the optical band. 

Since the effect of the speciflc window function is even 
stronger for X-ray data than for the optical data, we again 
performed Monte Carlo simulations, as described in Ap- 
pendix A. We considered two models for the shape of the 
NPSD. 

The first model was a broken power law identical to the 
adopted model for the optical power spectrum (see Eq. ISJ, 
with the frequency break and the slope above it being the 
free parameters of the model. The slope obtained from sim- 
ulations (see Table|5J is slightly flatter than the optical one, 
but the difference is within 1 a error. However, the best 
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Table 2. The NPSD parameters determined from Monte Carlo simulations for a single break model of 1910-2000 B light curve and two 
models of the X-ray light curve, one with a single and another with a double break 
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Quoted errors are determined from xlist statistics (see Appendix A) and correspond to 1 cr error for one parameter of interest. 



fit break frequency of the NPSD in the X-ray band corre- 
sponds to a timescale of around ^ 100 days. Formal contour 
errors slightly overlap (see Fig. 0J , particularly at the 90% 
confidence level but the probability maxima are well sepa- 
rated. The marginal probability plots made by integrating 
the probability distribution along one of the parameter axis 
(see the intermediate panel of Fig. |1J show that the distri- 
butions of the slopes in the optical and X-ray NPSD are 
only slightly different, namely they are shifted by ~ 0.4 
decades, while the distributions of the characteristic vari- 
ability timescale are essentially different. The value on the 
order of 100 days is favored for X-ray data but for the opti- 
cal data the maximum coincides with the longest timescale 
under study (~ 100 years). It strongly supports the view 
that the X-ray and optical variability properties at long 
timescales are significantly different. 

The second model was the shape frequently adopted 
for the galactic sources, a power law distribution with two 
breaks but fixed slopes of 0, 1 and 2 at the lowest, inter- 
mediate and the highest frequencies. Such a model is also 
a good representation of the observed power spectrum, al- 
though formally, the best fit model has a lower acceptance 
probability than in the previous case. The break frequencies 
in the best fit model correspond to the timescales of ~ 6 
days and ~ 1000 days. The low frequency break is at lower 
frequencies than expected from Fig. |S] since apparently the 
extra power at the lowest frequencies shifts it towards longer 
timescales. 

The spectrum may perhaps be better represented by 
a number of Lorentzian components, as suggested e.g. by 
Pottschmidt et al. (2002) for galactic sources but the quality 
of our light curves was not high enough to attempt more 
sophisticated fits. 



4.4 Structure function 

The analysis via the structure function was applied to the 
Crimean data covering the years from 1968 till to 2000 in 
UBVR bands, to the combined photographic plus photo- 
electric data in B band covering the years from 1910 till to 
2000, to the AGN Watch data, and also to the X-ray data 
collected from the literature. 



4^.4.1 Optical hand 

Fig.|S|shows SFs for the optical flux in various bands, where 
the Crimean data is plotted with filled circles. One can see 
that the shape of the SF which would have been inferred 
from the Crimean V magnitude data alone can be very well 
described as a single power law form, with the slope b « 
0.7 for time scale of variability from 1 to 3200 days. Such 
value of b shows that the dominating process is a single 
shot noise process. But in the U and B bands, the shape 
of the SF is more complex, and a clear break in the SF is 
seen at about 30 - 50 days, with much steeper slope at short 
timescales and a fiatter one in long timescales. The results of 
Merkulova et al. (2001) show that the steep slope continues 
down to the timescales of minutes. It can again indicate that 
in these two time intervals two different processes operate. 
We therefore measured separately the slope in time interval 
from 1 to 40 days and in time interval from 80 to 3200 
days (see Table O. Again, in the V band, a single slope of 
0.7 provides a good representation of the variability in all 
timescales shorter than 3200 days (8.8 years) but at shorter 
wavelengths the two slopes are significantly different. This 
difference is most visible in U band. The traces of similar 
structure are only barely seen in V band. 

We treated separately the SF of AGN Watch data. Con- 
tinuum fluxes at A 4600 A were converted to fluxes in B band 
and plotted as open circles in the third panel of Fig.jSJ. The 
corresponding points match well the Crimean data above the 
timescales of ~ 10 days but they steepen significantly be- 
low. This effect is probably due to the presence of a strong 
linear trend in the optical data during the period covered 
by AGN Watch monitoring (see Fig.0. Such a linear trend 
translates into SF oc r^, as can be shown analytically from 
Eq. 01 This kind of discrepancy underlies the importance of 
using as long time series of observations as possible in order 
to avoid artifacts. 

One of the important characteristics of the structure 
function is the point which marks the beginning of plateau at 
the long time lags. This time scale gives the maximum time 
scale Tmax of correlated signals or, equivalently, the mini- 
mum time scale of uncorrelated behavior. From Fig. |H] we 
can see that the SF reaches plateau at logr « 3.2 — 3.6, i.e. 
about 1600 - 4500 days. We can estimate it most accurately 
in the B band, due to the longest time coverage in this band. 



© 0000 RAS, MNRAS 000, 000-000 



Variability of accretion flow in the core of the Seyfert galaxy NGC 4151 9 



-1 



X-ray band 




-2 

log(f) [1/day] 




-3 -2.5 -2 
log fi [1/day] 



Figure 5. A plot of NPSD multiplied by the frequency, illus- 
trating the frequency regime where most of the power is present: 
dotted line - B band (1910 - 2000), continuous line - 2 - 10 keV 
band. Lower panel shows the error contours for the fit of the 
NPSD with two broken power laws of slopes 0, 1 and 2, and the 
break frequencies given by log /i and log /2 . Cross marks the best 
fit position. 
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Figure 6. The structure function for NGC 4151: upper panel - X- 
ray band, second panel - U band, third panel - B band (Crimean 
data - filled circles; AGN Watch - open circles; photographic data 
- triangles), lowest panel - V band. Errors given here are statistical 
errors (see Appendix B). Dashed line shows the level of SF due 
to the measurement error in the data. 



Since the photographic data is of relatively low accuracy, 
and the number of data points is small in comparison with 
photoelectric measurements, we compute separately the SF 
only for historical data. This structure function is shown in 
the third panel of Fig. |^ with open triangles. We see that 
the slopes of structure functions derived from both data sets 
are similar within the errors. Although the length of histor- 
ical data set is about 25200 days (logr = 4.4), the plateau 
begins earlier, at about log r — 3.2 (corresponding to the 
time interval 1600 days). So the existence of the plateau is 
not an artifact of an insufficiently long data train, but its 
presence is supported by long observations. 

The flattening observed in SF at 1600 - 4500 days (4 - 
12 years) corresponds well to the flattening hinted by the op- 
tical NPSD, at 10 years. The change of the slope at some 
temporal frequency is also seen in both cases, although it 



does not happen precisely at the same timescales. SF sug- 
gests a break at ^ 30 while that apparent in the NPSD is 
at ~ 500 days. 

Although the SF is less sensitive to the window func- 
tion problem than the NPSD technique, we performed some 
Monte Carlo simulations to estimate better the uncertain- 
ties involved in the analysis. We adopted a shot noise model 
of the light curve, with the model parameters 

dT, a, P, b^2 + a + 2P, (10) 

being the mean time separation of the flare, power law in- 
dex of flare number distribution, power law index of flare 
duration distribution, minimum and maximum duration of 
a flare, correspondingly (see Appendix B). 

The best fit model parameters representing the Crimean 
data in the U, B and V bands are given in Table and 
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Figure 7. Upper panel: the structure function in three color 
bands for the Crimean data (U - continuous line, B - thick dashed 
line, V - thin dashed line), with the contribution from observa- 
tional errors subtracted. Lower panel; comparison of the SF for 
the data (dashed line) and for the best fit model (thick line) in 
the V band. The error on the simulated SF in the lower panel is 
determined from the dispersion in Monte Carlo simulations (see 
Appendix B). Dotted lines show the level of SF due to the mea- 
surement error in the data. 

the best fit models are shown in Fig. [7| There is a good 
overall agreement between the SF for the best fit model 
and the SF derived from the data (see lower panel of Fig. Q 
for comparison of the observed and modeled SF in V band). 
The SFs in all three colors are basically similar to each other 
although hints of systematic differences may be seen at the 
intermediate timescales between 1 and 60 days, reflecting the 
minor difference in the best fit parameters. Those differences 
are marginally significant. In the lower panel of Fig. Q we 
show the dispersion of the SF for the family of Monte Carlo 
curves in the V band; the indicated error is about 0.2 in 
logarithmic scale at intermediate timescales and becomes 
larger at long and short timescales. 

We also used all the data collected in the B band in 
order to exploit fully the data length. The result is shown 
in Fig. |S] In this Figure we plot the normalized structure 
function for better comparison between various curves, and 
we also subtract the value 2a1,.^, i.e. twice the variance of the 
measurement uncertainties. On the old photographic data 
we adopt a measurement error of 8.5%. We see that the SF 
is not significantly different from the one obtained from the 
Crimean data alone. However, the apparent flattening of the 
SF at the longest timescales is still better visible due to the 
access of the longest timescales. 

An interesting question is whether the variability prop- 
erties during flares are different, and to address this, we 
repeat our analysis for the two cycles A and B separately. 
All slopes of the SF are presented in Table |21 An inspec- 
tion of TableOreveals that the properties of the light curves 
{SF slopes) are significantly different in the two cycles. Dur- 
ing the cycle B the slope difference between the short and 
long timescales is moderate, although a feature around ~ 50 
days is present. The plateau for long time intervals begins 
at logr=3.2 - 3.4. During the cycle A, the difference of the 
slopes is apparent in all color bands, with the long timescale 
slope being very flat (Fig. |U] see also Fig|H] second panel). 



-2-101234 




-2-10 1 2 3 



log T (days) 

Figure 8. Normalized SF for observational data. Upper panel: 
results for X-ray band with long timescale trends determined from 
daily averages, and short timescale trends from the ASCA data 
alone, together with the SF for the cycle A averaged over UBV 
band. Second panel: comparison of SF averaged over UBV bands 
for the cycles A and B. Third panel: comparison of the result 
from all combined photographic (ptg) and photoelectric (phe) 
data in B band with photoelectric Crimean B band data alone. 
Lowest panel: comparison of the X-ray SF with V band SF based 
on Crimean data (see Fig. |7J and best fit V band model; the 
envelope of the V curve marks the error estimate from Monte 
Carlo simulations. 



suggesting that long timescale trends are relatively weak. 
Monte Carlo simulations with short lasting shots (20 days 
and less) reproduced well the observed light curve proper- 
ties. This is caused by the fact that in the cycle A we do not 
see a single large outburst but instead, a number of weaker 
eruptions superimposed on a weak longer trend. The local- 
ization of Tmax, which is the time scale when SF reaches 
plateau, is less clear in the cycle A than in the cycle B, or in 
the composite data. During the cycle A there are no signif- 
icant differences in the SF shape for all three UBV colors. 

The systematic difference between the cycles is also re- 
flected in our Monte Carlo simulations (see Table^J- We see 
most clearly a shift (reduction) in the maximum duration 
of the flares in the best flt model for the cycle A alone in 
comparison to the model for the whole Crimean data: this 
reduction is from 800 days to 20 days. 
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Table 3. Structure function slopes b in various color bands, data sets, and timescales. 



Band 


b 




r 


b 


a{b) 


r 



Cycle A (25.03.1968 - 03.07.1988) 

(1 - 40) days (80 - 3200) days 



u 


1.24 


0.08 


0.95 


0.15 


0.04 


0.61 


B 


1.15 


0.10 


0.94 


0.18 


0.04 


0.66 


V 


1.07 


0.10 


0.93 


0.17 


0.04 


0.64 


mean UBV 


1.04 


0.08 


0.95 


0.16 


0.03 


0.64 



Cycle B (11.02.1989 - 13.03.2000) 

(1 - 40) days (80 - 1500) days 



u 


1.24 


0.05 


0.98 


0.58 


0.04 


0.94 


B 


1.12 


0.03 


0.99 


0.64 


0.04 


0.96 


V 


0.86 


0.04 


0.98 


0.67 


0.04 


0.96 


mean UBV 


1.02 


0.04 


0.99 


0.63 


0.04 


0.96 



Cycles A + B (25.03.1968 - 13.03.2000) 

(1 - 50) days (80 - 3200) days 
U 1.21 0.08 0.96 0. 60 0.03 0.97 
B 1.06 0.06 0.97 0. 64 0.03 0.98 
V 0.95 0.05 0.97 0. 66 0.02 0.98 



Photographic observations + Cycles A + B (5.03.1910 - 13.03.2000) 

(1 - 40) days (80 - 3200) days 

B 0.99 0.07 0.95 0. 60 0.03 0.98 



X-rays (20.10.1974 - 05.01.2000) 

80 min - 10 days (10 - 4700) days 

2-10 keV 1.24 0.04 0.98 0. 23 0.02 0.90 

Spectral bands are marked in column (1), time interval for correlated variability are specified above columns (2 - 4) and above columns 
(5-7), slopes b and errors of slope, iT(b), are shown in column (2), (3) and (5), (6), and correlation coefficients r for the corresponding 
linear regressions are in column (4) and (7). Errors of the slopes are determined from fitting the power law to SF points in a specified 
timescale range. 

seen at logr around -1.8 and 0, and a complex structure is 
visible at timescales of 10 days. 

Because the scatter was so large, we reanalyzed the X- 
ray data as follows. We first created daily averages from all 
X-ray observations and computed the long timescale part 
of the SF. Next we computed the short timescale SF from 
the ASCA long monitoring alone, up to logr = 0.8 (i.e. 
timescale corresponding to a half of the total duration of the 
ASCA 2000 monitoring). We renormalized the SF dividing 
it by the variance, (see Eq. (HJ. Both NSFs matched each 
other well at the transition timescale. The results are shown 
m Fig. IHl The new plot appears smoother, but the overall 
shape did not change. 

X-ray SF shows clearly the change of the slope around 
the timescale of a few days and a trace of flattening at the 
longest timescales. There is therefore a good correspondence 
between the SF and the NPSD in the X-ray band. There is, 
however, a systematic difference between the SF in the X-ray 
band and in the optical band, hinted in our NPSD analysis. 
The difference is seen most clearly when all optical data are 
included (see the lowest panel of Fig.|H|for a comparison with 
V band). However, the difference is much smaller when only 
the photoelectric measurements from cycle A are taken (see 
the second panel of Fig.jHJ. In this case a significant differ- 
ence is seen only at intermediate timescales from a fraction 
of a day to 10 days. The presence of the two characteristic 
timescales in X-ray band - ~ 500 days and ~ 5 days - should 
be confirmed by continued monitoring of this source. 



Table 4. The SF model parameters determined from Monte Carlo 
simulations. 



Cycle A A + B 



Band U U B V 



dT [days] 0.10 0.20 0.25 0.25 

13 0.45 0.25 0.25 0.17 

a -1.47 -1.80 -1.70 -1.60 

Tmin [days] 0.01 0.01 0.01 0.01 

Tmax [days] 20 800 800 800 

log(T)„i„ -0.05 0.03 -0.05 0.03 

\oglT)„^ax 1.3 3.0 3.0 3.0 

b 1.22 0.70 0.79 0.74 

error b (Icr) ±0.17 ±0.07 ±0.06 ±0.06 



The quoted values dT, a, (3, Tmin a^nd Tmax are best fit model 
parameters for flare separation, power law indices in the number 
of flares and flare duration distributions, and the minimum and 
the maximum flare duration. The slope b of the SF was next de- 
termined in the time interval between log(r)min and log(r)maa;- 



4-4-2 X-ray band 

The shape of the SF in the X-ray band can be approximately 
represented by a single power law with a slope ~ 0.3 (see 
Fig-inj- A fiattening at the timescales above ~ 1000 days is 
not excluded but it is not strongly required. Errors in this 
part of the plot are large. Minor changes of the slope are also 
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Figure 9. Optical light curve properties during the cycle A. Up- 
per panel: the observed light curve in U band and a representative 
light curve from Monte Carlo simulations, normalized to the mean 
flux and offset for clarity. Middle panel: observed SF in three color 
bands. Lower panel: mean and dispersion of the SF from Monte 
Carlo simulations for best fit model parameters, together with the 
observed SF in the U band. 



4.5 Autocorrelation function 

We supplement the analysis by calculating the autocorre- 
lation function for the full data set in the U band and in 
the X-ray band. The difference between the two plots is re- 
markable. We see from Fig. 1101 that X-ray flux shows very 
good correlation only at extremely short timescales below 10 
days. At timescales of order of 10-1000 days, corresponding 
to a hump in a X-ray NPSD (see Fig. |^ , a broad shoulder 
replaces the sharp peak and falls down slowly. A significant 
minimum in the correlation is seen at a timescale of ^ 215 
days, which, when multiplied by 2, corresponds well to the 
first minimum in NPSD. It is consistent with the SF not be- 
ing completely flat at such a long timescales (see Table |3J|. 
The correlation reaches zero at 400 days, giving timescales 
of 800 days, in agreement with the presence of considerable 
power in NPSD down to timescales of ~ 1000 days. 

U band autocorrelation function does not show a similar 
structure and falls down steadily, crossing at a time lag 
of 1650 days and indicating the characteristic timescale of 
about 3300 days. This is again in rough agreement with the 
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Figure 10. Autocorrelation function in U band (continuous line) 
and in the X-ray band (dashed line). 

determination from SF (- 3200 days) and NPSD (~ 5000 
days). 



5 DISCUSSION 

5.1 Character of optical and X-ray variability 

Our results indicate that there are two variability mecha- 
nisms operating in the Seyfert galaxy NGC 4151, and we 
infer this by comparing the results of our analysis of the 
optical and X-ray data. 

The NPSD plots for the optical band suggest the flat- 
tening at the timescales of ~ 10 years and Monte Carlo sim- 
ulations show that this value is actually a lower limit for the 
characteristic timescale. The SF plots for the optical band 
suggest a flattening at ~ 4 — 10 years and the Monte Carlo 
simulations show that models with shot duration shorter 
than 800 days cannot explain the data. 

On the other hand, the NPSD plots for the X-ray band 
suggest that most of the power is present in the intermedi- 
ate timescale range, between 5 and 1000 days. Monte Carlo 
simulations support this view, giving the best fit flattening 
timescale at 130 days and a 1 cr error upper limit of 5 years. 
SF plots in the X-ray band show that the power rises mostly 
at 1 - 10 day timescale range. The distribution is rather fiat 
at longer timescales, with a hint of flattening above 1000 - 
2000 days. 

In summary, although errors are large, we conclude that 
the optical variability occurs predominantly on timescales 
longer than ~ 10 years while the X-ray variability is mostly 
confined to the timescale range 5 days — 3 years. Therefore, 
two different variability mechanisms are apparently at work. 
Our results therefore support the early claim of Belokon, 
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Babadzhanjantz, & Lyutyi (1978) based on much more lim- 
ited data. 

The difference in the dominating timescales in the op- 
tical and X-ray band does not exclude some interaction be- 
tween the two bands. Long time scale behavior of the vari- 
ability of X-ray flux appears to be somewhat correlated with 
the optical flux, suggesting an influence of the structure of 
the optical emitter on the X-ray emitter. Also, the fast vari- 
ability, similar to the X-ray variability, is present in the op- 
tical data which possibly reflects the effect of the X-ray ir- 
radiation of the cool material. 

The shortest characteristic timescale seen in our X-ray 
data is in agreement with the timescale 13^5^ days deter- 
mined by Collier and Peterson (2001). We do not see clearly 
any change of the NPSD slope if we try to extend our com- 
putations towards shorter timescales but the SF shows a 
change of the slope just above this range. Collier and Pe- 
terson (2001) did not flnd the long timescale component be- 
cause their data for NGC 4151 covered only 99 days. Their 
conclusion about the overall similarity between the opti- 
cal/UV variability and X-ray variability on short timescales 
agrees with our results, although our results are not partic- 
ularly accurate in that range. 

Our inference that the long timescales play the domi- 
nant role of the optical band cannot be an artifact of the 
complex analysis since it is actually apparent from Fig. 
The actual value of the characteristic timescale is strongly 
influenced by the recent outburst (cycle B) lasting from 1989 
till 2000, although hints for similar timescales are present in 
the older data as well, when they are superimposed on low 
amplitude longer trend. 

The lack of trends longer than ~ 3 years in the X- 
ray data obtained from NPSD is less secure and a clear ab- 
sence of longer timescales is also not immediately apparent 
from the SF analysis. Large errors associated with the early 
measurements make this determination rather difficult. A 
few more years of further monitoring will resolve this is- 
sue definitively. However, the dominant role of much shorter 
timescales in the X-ray band seems to be well established. 



5.2 Correlation of optical and X-ray luminosity 

Long data sets provided also an opportunity to better study 
the correlation between the U band and the 2-10 keV X- 
ray emission. The correlation between the whole data sets 
is surprisingly poor. Short timescale study of UV and X-ray 
variability of Edelson et al.(1996) showed a very strong cou- 
pling during 10 days of monitoring. Noticeable correlation 
was found by Perotti et al. (1990) between the hard X-ray 
35 keV flux and the optical flux. 

Oknyanskij (1994) discussed the issue of the dependence 
of broad band correlations on the luminosity state of the 
nucleus and he found that in the low state in 1983 - 1986, 
correlations are more significant than in the high state in 
1975 - 1978. The same conclusion was drawn by Perola & 
Piro (1994). 

X-ray and U band measurements in our data sets are 
frequently significantly separated in time and the old X- 
ray data are not very accurate. Therefore, we repeated our 
analysis for a subset of data points including only the rela- 
tively new and accurate measurements (observations start- 
ing from 1990) and taking only those pairs for which the 
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Figure 11. The correlation between the optical flux in U band 
and the X-ray flux for the data points collected starting from 
1990, with U and X measurement separation smaller than 1 day. 



separation of the measurement time is smaller than 1 day. 
Our data sets contain 37 such pairs. The result is shown 
in Fig. 111! The distribution observed in the plot is consis- 
tent with the coefficient of correlation r = 0.780, and is 
highly inconsistent with the null hypothesis Ho assuming 
no correlation, i.e. r — 0. We tested it following the pro- 
cedure of Borczyk,Schwarzenberg-Czerny & Szkody (2003). 
For n = 35 degrees of freedom, our value of r tested against 
Ho yields Student t = 7.36 corresponding to the significance 
level 1 - 2.6 • 10"^ 

The result strongly confirms the trend seen by Oknyan- 
skij (1994) and Perola & Piro (1994): the correlation is 
stronger when the source is fainter, and the X-ray flux shows 
significant systematic 'deficit' when the source is bright. It 
shows that the smaller X-ray rms variability (see Table 
does not result from the presence of a constant X-ray com- 
ponent. In opposite, we rather have a small constant U com- 
ponent (offset), a linear trend for moderate fiux and a flat- 
tening for high flux, but with significant dispersion. 



5.3 Comparison with NGC 5548 

Although the power spectra in the X-ray band are available 
now for a number of AGN, the variability in the optical 
band was quantitatively studied only for the Seyfert galaxy 
NGC 5548. Therefore, we can only compare the results for 
the two Seyfert galaxies: NGC 4151 and NGC 5548. 

Below, we compare the structure functions in U, B, V, 
Rj bands of NGC 4151 with the resufis for NGC 5548 ob- 
tained by Doroshenko et al. (2001). Our result is that the SF 
of NGC 4151 in U and B bands has a break (see Table|21and 
Figures |S| , whereas SF of NGC 5548 has a single slope 
in all bands. The values of the slope, 6, are 0.72 ± 0.02 
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for U band and 0.69 ± 0.02 for B band roughly similar to 
our results for the long timescale behavior in Cycles A-'rB. 
However, the timescale where the SF flattens may be some- 
what different in the two objects: in NGC 4f 5f it is about 
logr « 3.2 — 3.6 [day] while in NGC 5548 it is somewhat 
shorter, just below logr = 3.0 [day]. Note, however, that the 
slopes for NGC 5548 were determined without correcting for 
the level of the observational error. 

We can also compare the normalized power spectrum 
density of NGC 4151 in the optical band (Fig. ^ with 
NGC 5548 given by Czerny et al. (1999) in their Fig. 7. 
We can see an apparent maximum in Power x Frequency 
representation in both objects, at log / ~ —3.7 [1/day] for 
NGC 4151 and at log / ~ -3 for NGC 5548. However, Monte 
Carlo simulations show that for NGC 4151 we actually see a 
lower limit of the dominant timescale and the same can be 
true for NGC 5548 with data set covering shorter period. 

We also compare the normalized power spectrum den- 
sity of NGC 4151 in the X-ray (2 - 10 keV) band (Fig.lHJ with 
NGC 5548 given by Uttley et al. (2002). The overall shape 
of PSD in both objects is similar. The PSD of NGC 5548 is 
well represented by a single power law (q = 1.6) up to log 
/ = -3. For NGC 4151 the fit is better for a solution with a 
break (a = 1.5, log / = -2.1) but a single power law is also 
acceptable. 

The masses of the central black holes estimated from 
reverberation studies are log M = 7.83 for NGC 5548 (Pe- 
terson & Wandel 1999) and logM = 7.08 for NGC 4151 
(Wandel, Peterson, & Malkan 1999). Mass determination 
from accretion disk fitting gives log M = 7.4 for NGC 5548 
and log M = 7.0 - 7.3 for NGC 4151 (Czerny et al. 2001). 
The high frequency tail of our NGC 4151 power spectrum 
is shifted down by a factor of 1.5 in comparison with the 
NGC 5548 power spectrum of Uttley et al. (2002). Scaling 
used by Hayashida et al. (1998) in such a case implies that 
the black hole mass in NGC 4151 is higher, in contrast to 
the above results. This suggests that indeed the black hole 
masses in NGC 4151 and in NGC 5548 are comparable, but 
again, the errors are large, probably higher than quoted in 
the above papers. 

5.4 Comparison with Cyg X-1 

The comparison of observational characteristics of galactic 
black hole candidates (GBH) and AGN is not straightfor- 
ward. Observed photon spectra generally consist of two prin- 
cipal spectral components: soft bump usually interpreted as 
the accretion disk emission, and a power law interpreted as 
a result of Comptonization by a hot optically thin plasma. 
In AGN these components are well separated: the first one 
dominates the optical/UV emission, while the second one 
dominates the X-ray band. In galactic sources both are seen 
in X-ray band but their dominant role depends on the lu- 
minosity state of the source. Therefore, the power spectrum 
of a GBH in the hard, power-law dominated state should be 
compared with the X-ray power spectrum of an AGN while 
the power spectrum of a GBH in the soft disk-dominated 
state should be compared with the optical power spectrum of 
an AGN. While the former comparison has been attempted 
often in the past, the latter has not been addressed exten- 
sively. 

X-ray power spectra of AGN are generally quite simi- 




Figure 12. The NPSD for Cyg X-1 in the soft state, June 1996, 
(filled squares, from Gilfanov et al. 2000) and in the hard state 
February 3, 1997, (open circles), together with the NPSD for 
NGC 4151 in the optical band (dotted histogram) and in the 
X-ray band (continuous histogram), both shifted horizontally by 
6.3 to account for a mass difference between the two sources. 



lar to power spectra of galactic objects in their hard state 
(see e.g. Czerny et al. 2001 and the references therein). The 
power peak in / x NPSD diagram is broad, covering about 
2 decades, and its position on the frequency axis depends 
linearly on the mass of the black hole and can well be used 
for the mass determination (Hayashida et al. 1998; Czerny 
et al. 2001). A comparison of the X-ray power spectra of Cyg 
X-1 and NGC 4151 (Fig. I12I I generally support this picture, 
and suggest that the mechanisms responsible for the emis- 
sion and variability of the power law spectral component are 
indeed identical for AGN and GBH. 

However, a comparison of the optical power spectrum of 
an AGN with power spectra of galactic sources in soft state 
does not show any striking similarity. In galactic sources 
the power is evenly distributed across the broad range of 
frequencies. In AGN, the longest timescales dominate. Our 
results for NGC 4151 confirm such a view (see Fig. [T^ . 
Timescales dominating the optical/UV variability in AGN 
(years) correspond to the timescales of hundreds/thousands 
of seconds in galactic sources. Strong outbursts with such 
timescales are only seen in the microquasar GRS19154-105 
(see Belloni et al. 1997, Janiuk, Czerny, & Siemiginowska 
2000). However, this source shows strong jet outflow and 
therefore it is not a simple analogue of a radio quiet Seyfert 
galaxy. 



5.5 Decomposition of optical light curve into fast 
and slow variability components 

The presence of the two separate types of variability - a slow 
component and a fast component - was already suggested by 
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Belokon et al. (1978). Our results from the analysis of the 
PSD and SF in U and B bands based on much better data 
coverage confirm this result. Therefore, the data are not well 
represented by a sum of a constant component and a rapidly 
variable component, as suggested by Ulrich (2000). 

In order to visualize the presence of the two sepa- 
rate types of variability - a slow component and a fast 
component- we complement the computations of the optical 
NPSD and the SF with a simple illustration. We calculate 
the smooth long timescale light curve at each data point by 
averaging all available data points separated from a given 
instance by less than adopted smoothing timescale, AT. A 
representative set of such light curves for a range of AT is 
shown in Fig. 

Adopting AT — 500, we display the short timescale 
variations by subtracting the smoothed light curve from the 
original one (see Fig. 1181 panel (a)). The result is shown in 
Fig. 1131 panel (b) . We see that the amplitude of variations 
display significant trends. This dependence of the amplitude 
on the fiux was discussed by Doroshenko et al. (2001). 

Such trends almost disappear if we plot the normalized 
light curve, i.e. if the value of the flux difference is divided 
by the value of the flux of the smoothed light curve (see 
Fig. 1131 panel (c)). This figure provides an illustration of the 
same fact that was stressed by Uttley & McHardy (2001): 
only the normalized power spectrum is independent from 
the luminosity state of the source. 

We can now compare these trends with variability of 
Cyg X-1 in its soft (disk-dominated) spectral state. In 
Fig. 1131 (panel (d)) we show a sequence from the light curve 
of Cyg X-1 in February 3, 1997 (30 yr for NGC 4151 are 
equivalent to 300 s of Cyg X-1). We see that this light curve 
is clearly different from the original U-band light curve of 
NGC 4151 (see Fig. I^J but there is some similarity of Cyg 
X-1 light curve to that shown in panel (c) where all the long 
trends from NGC 4151 light curve have been removed. This 
supports our earlier conclusion that the fast optical vari- 
ability component in AGN may well be connected with the 
X-ray reprocessing, and the X-ray variability itself is similar 
in AGN and GBH so the nature of fast variability might be 
the same in both types of objects. 

5.6 Intrinsic variability versus obscuration 

A number of authors suggested that the variable obscuration 
is, at least partially, responsible for the observed variability 
of AGN (e.g. Collin et al. 1996, Boiler et al. 1997, Brandt et 
al. 1999, Abrassart & Czerny 2000, Risaliti, Elvis, & Nicas- 
tro 2002). 

The absorbing column density in Ginga/EXOSAT data 
seems to be lower when the source is brighter (Yaqoob & 
Warwick 1991). However, if we look at the most recent and 
accurate observations, no obvious correlation is seen. In De- 
cember 1993 when the AGN Watch monitoring was done the 
source was bright in U band ('^ 75 mjy; Kaspi et al. 1996, 
Crimean data) and the intrinsic X-ray column was 3.5 x 10^^ 
cm"^ (Edelson et al. 1996). In the beginning of March 2000 
the source was faint in the U band (~ 21.5 mJy; Crimean 
data) and the intrinsic X-ray column determined from the 
Chandra data was 3.7 x 10^^ cm"^ (Ogle et al. 2000). 

On the other hand, the observed delays between the 
optical continuum and the emission lines (e.g. Kaspi et al. 



1996) and between the optical emission and the infrared 
emission (e.g. Oknyanskij et al. 1999) show that the vari- 
ability on the timescales of days is certainly intrinsic to the 
source, at least partially. 

lUE monitoring during 18 years, which covered both the 
very low state around 1988 and the luminosity peak in the 
middle of outburst B around 1995 showed a factor over 12 
increase in the flux although the UV line CIV1550 showed 
lower amplitude, by a factor of 5 (see Fig. 11 of Ulrich 2000). 
It might suggest that the continuum variations are addition- 
ally enhanced by some variations in the absorbing medium, 
perhaps due to variable ionization. 

This shows that although the variable obscuration can- 
not be totally excluded, it cannot be the only factor respon- 
sible for the observed variability and most of variations are 
instead intrinsic to the source. 

5.7 Interpretation of the variability 

The strong variability of AGN in optical, UV and X-ray 
band is a well known but not well understood phenomenon 
(see Mushotzky, Done, & Pounds 1993; Ulrich, Maraschi, & 
Urry 1997). The observed properties of the X-ray time series 
indicate some stochastic process (McHardy & Czerny 1987; 
Czerny & Lehto 1997), which physically can be interpreted 
as magnetic flares above the cold disk or shocks developing in 
the hot accreting material. Optical variations are partially 
caused by reprocessing (although modeling is not simple; 
see e.g. Rokaki, Collin-Souffrin, & Magnan 1993; Kazanas 
& Nayakshin 2001; Zycki & Rozahska 2001; Wang, Wang, & 
Zhou 2001) and partially by another process operating on 
longer timescale, possibly the radiation pressure instability 
(e.g. Czerny et al. 1999). 

The variability of the flux from the cold disk can be 
related to the Keplerian timescale, tx, thermal timescale, 
tth or viscous timescale tvisc of the accretion flow. The first 
value is universal (and depends only on the mass of the 
central object), and the other two can be computed for 
a Shakura-Sunyaev disk model. The variations due to the 
changes of properties of the hot plasma may also be related 
to tK and tth, although the lack of theory for the dynamics 
of the hot phase makes this parameterization less firm. 

To compute these time scales, we fix the mass of the 
black hole at 3.7 x 10^ Mq, based on variability studies, and 
take the average luminosity to the Eddington luminosity ra- 
tio 0.02 resulting from the black hole mass value and the 
estimate of the total bolometric luminosity (see Czerny et 
al. 2001). For those parameters, we obtain the radial depen- 
dences of those timescales: 

tA' = 0.18m3.7ri/'' [d], (11) 
tth = l.Sao.lma.rrli" [d], (12) 

tmsc = 4000Qo.im3.7ri^^m^02 [d], (13) 

where r is expressed in units of 10 Schwarzschild radii, m in 
units of 3.7 x IO^Mq, and a standard viscosity parameter 
a in units of 0.1. The mass of the black hole enters the for- 
mulae linearly, and the accretion rate affects only the third 
quantity. In these computations we adopted MEdd = 1-69 x 
10^** X {M/Mq) g s-\ and Ledd = 1-27 x 10^^ x {M/Mq) 
erg s~^, corresponding to the Newtonian accretion efficiency 
of 1/12. 
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Figure 13. Panel (a): the smoothed light curve of NGC 4151 in U band, where smoothing timescales equal 200, 500 and 1000 days were 
used. Panel (b): the light curve of NGC 4151 in U band, with the long timescale trend of 500 days subtracted. Panel (c): the normalized 
light curve of NGC 4151 in U band, with the long timescale trend of 500 days subtracted. Panel (d): X-ray light curve of Cyg X-1 when 
the source was in the soft state (May 30, 1996). Errors in the Crimean data are ~ 2.2 %, errors in Cyg X-1 light curve arc ~ 25 cts . 



We can now compare those timescales with the char- 
acteristic points present in the NPSD and SF both in the 
optical and in the X-ray band. The derived timescales are 
related to the Fourier frequencies: / = (27rt)~^. The fac- 
tor 2-K is important if we want to compare the predicted 
timescales with the characteristic frequencies in the shape 
of the NPSD. Therefore, observed characteristic frequencies 
(1/5, 1/1000 and a lower limit of 1/4000 d"^) translate into 
0.8 d, 160 d, and more than 640 d, respectively. Here we used 
the timescales resulting from two-break fit (see Tablel^J, mo- 
tivated by the analogy with the galactic sources. 

The shortest X-ray timescale corresponds to the ther- 
mal timescale at the marginally stable orbit of the non- 
rotating black hole if the viscosity parameter is ~ 0.03, 



as obtained from spectral fits to the data for NGC 5548 
(Kuraszkiewicz, Loska, & Czerny 1997). The value of the 
longest X-ray timescale is by a factor of 200 longer which 
suggests that the X-ray production region extends in that 
case up to ~ IQQRschw , if timescale scaling as r^^'^ holds. 

The timescale dominating the optical variability seems 
to be too long for pure thermal variations related to the 
radiation pressure instability. The thermal timescale of 640 
d corresponds to a distance of ~ 220 Rschw (assuming a 
as above), much larger than the extension of the instability 
zone (~ 70 Rschw, calculated using the code of Roanska 
et al. 1999) for the adopted system parameters. Viscous 
timescales (see Eq. I13II on the other hand, are much longer 
than the period covered by the optical data. However, re- 
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cent computations performed by Szuszkiewicz (1999) indi- 
cated much faster evolution than estimated from Eg. 1131 - an 
outburst lasted ~ 11 years for adopted model parameters: 
W^Mq black hole, a — 0.1 and accretion rate m — 0.06. 

Therefore, the observed variability in the optical band 
neither obviously proves nor contradicts the supposition that 
the optical variability at long timescales is caused by the ra- 
diation pressure instability of the optically thick disk. Fur- 
ther optical monitoring as well as theoretical progress are 
needed to resolve this issue. However, the lack of strong 
long timescale (hundreds of seconds) variability and the 
lack of radiation pressure dominated zone in Cyg X-1 (see 
e.g. Gierlihski et al. 1999; Janiuk, Czerny, & Siemiginowska 
2002), when compared with the presence of long timescale 
variations and the disk unstable zone in NGC 4151, favors 
the radiation pressure scenario. Analysis of the optical vari- 
ability of NGC 5548 supported this view as well (Czerny et 
al. 1999). 



6 CONCLUSIONS 

Our analysis of the 90 years of the optical data and 27 years 
of the X-ray data for NGC 4151 using the NPSD and the 
SF techniques, gives the following results: 

• variability properties in the optical and X-ray band are 
significantly different 

• X-ray variations are predominantly in the frequency 
range 1/5 - 1/1000 

• optical variations are dominated by a component with 
a frequency on the order of (or possible smaller than) 
~ 1/5000 d~^ (or ~ 10 years), but the short timescale vari- 
ability may be related to X-ray variability 

• the presence of long timescale variability in NGC 4151 
and the absence of analogous variability (on timescales of 
hundreds of seconds) in Cyg X-1 favors the radiation pres- 
sure mechanism 

• the radial extension of the X-ray production region is 
similar in NGC 4151 and Cyg X-1 so the possible presence 
of radiation pressure dominated region in NGC 4151 and its 
absence in Cyg X-1 does not affect the X-ray production 
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APPENDIX A: DETERMINATION OF THE 
ERROR ON THE POWER SPECTRUM 

Two methods can be used to estimate the errors on the 
power spectrum determined from observations. The first 
method is to assign directly the errors on the plotted his- 
togram, and the second one is to do so through the Monte 
Carlo simulations. In this paper, we use both methods. 



A.l The direct method 

The direct measurement error is estimated as follows. Dis- 
crete power spectrum computed at Fast Fourier Transform 
(FFT) frequency grid has a X^(2) distribution. Its variance 
matches the expected power. Whether the same is true for 
unevenly sampled data remains yet to be proven (e.g. Tim- 
mer & Konig 1995). Below, we adopt the same procedure, 
formal objections notwithstanding. Since we use logarithmic 
values (see Papadakis & Lawrence 1993), we estimate the 
error on the logarithmic value as follows. If the true value 
of the power is distributed around the measured value uni- 
formly between and twice the value, one sigma negative 
error would correspond to a value equal to 0.16 of the mea- 
sured value. Therefore, the error of the logarithmic value 
on the negative side is equal to 0.43. We assume now that 
the errors are symmetric and adopt this value as a single 
measurement error of a single power spectrum in a single 
bin. 

The power spectrum is computed in the frequency range 
determined by the data. It is subsequently binned into log- 
arithmic bins of the width of 0.2. Measurements grupped 
into a single bin usually show significant dispersion, some- 
times even higher than the single measurement error as- 
signed above. Therefore we compute the error of a power 
spectrum in a bin by combining single measurement errors 
with the dispersion around the mean value x 

<^b{fk) = (2^ + ■ (14) 

When separate sequences are computed overlapping in 
the frequency range we add the spectra in the common fre- 
quency bins assuming that the weight of the spectrum is 
inversely proportional to the measurement error of the con- 
tributing spectrum. The error of the resulting spectrum is 
again computed by combining the errors of the single mea- 
surements al (determined as described above) with the dis- 
persion error 

'^hUk) = (2^ -2 + 2^ —2-) . (15) 

This method is very simple and does not depend on 
the model of the power spectrum. However, it does not in- 
corporate well the systematic errors caused by the window 
function. 



A. 2 Monte Carlo simulations of light curves 

This method allows to estimate also the systematic effects 
of the window function as well as to draw the error contours 
for the adopted parameterization of the shape of the power 
spectrum. 

In our approach we follow the method used by Uttley 
et al. (2002) and based on the method of Timmer & Konig 
(1995). It is equivalent to the method used by us in our pre- 
vious paper (Czerny et al. 1999) which consisted of filtering 
the white noise with the filter corresponding to the adopted 
power spectrum distribution. 

In our specific case we proceed as follows: 
• We assume a shape of the NPSD in a form of either a 
broken power law with the slope below the frequency break 
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(parameters: frequency break and the slope above it) or in 
the form of a power law with two breaks and the slopes 0, 1 
and 2 (parameters: two frequency breaks). 

• We generate a random realization of this distribution by 
assuming the random uniform distribution of the power 
spectrum around the adopted value and the random distri- 
bution of Fourier phases. We generate an artificial equally 
spaced light curve that is somewhat longer than the dura- 
tion of the data set through the FFT method. Next the light 
curve is folded with the observational window function thus 
reducing the long light curve to the number of points and 
spacing as in the observational data set. 

• Next, we follow exactly the procedure applied to the obser- 
vational data, i.e. we form yearly averages and other data 
subsets (see Sect. 13.111 . We compute the power spectrum 
NPSDi{f) as we did for the observational data. We repeat 
the entire procedure i = 100 — 500 times in order to create 
good statistics. The method is computationally intensive, so 
we cannot significantly increase this number for such a long 
and densely sampled data sets eis used in our paper. From 
this distribution, we determine the average simulated power 
spectrum NPSD{f), and the dispersion a{NPSD{f)). We 
can assign a Xdist(*) value to each of the simulated light 
curves: with the partial results forming the statistics 
around it 



2 



[NPSD'{fk)~NPSD{fk)Y 
a{NPSD{fk)r 



(16) 



The normalization of the power spectrum is adjusted to min- 
imize the Xdistii) value. 

• We also compute the same Xdist(o'^'S) value for the mean 
simulated power spectrum and the power spectrum deter- 
mined from observations: 



2 

Xdist 



_ [NPSD°'"'ifk) - NPSDjf, 



<j{NPSD{UW 



(17) 



• The quality of the fit for a given power spectrum pa- 
rameters is determined from the comparison of the value 
Xdisf(o6s) with the distribution given by Eq. 1161 the per- 
centage of the light curves with the Xdist(^) smaller than 
Xdist(o6s) gives the probability P of the model rejection. 

• Repeating the analysis for the range of the power spectrum 
parameters we find the best fit value (the highest acceptance 
probability) and we check if it is an acceptable model. 

• Next, we determine the error contours for the two parame- 
ters of interest. For that purpose, we assume that the proba- 
bility of any sets of values is proportional to the acceptance 
probability obtained above. We normalize the probability 
through performing the two-dimensional integral. We obtain 
the contour errors corresponding to the rejection probability 
of 68% (~ 1(7 error) and 90% by integrating this renormal- 
ized probability over the surface through subsequent prob- 
ability contours until the value of the integral is equal to 
0.68 or 0.9, correspondingly. Since the computations were 
performed on a regular grid, in practice we find for example 
the contour level Po.68 on the parameters of interest x and 
y iteratively from the expression 



0.68 = 



O .2 




(18) 



-3 -2 -1 

log(f) [1/day] 



Figure 14. Observed NPSD in the X-ray band (continuous his- 
togram) with errors determined directly from the data plotted 
against the best fit model (dotted histogram) of a single break 
power law with log/o = —2.1, a = 1.5. The xli^t for this 
model is 13.39 for 22 dof (25 frequency points, 3 model parameters 
including normalization). 



Such contours are slightly smaller than contours plotted di- 
rectly from the acceptance probability. 

The comparison of the observed NPSD with the best fit 
model in X-ray band is shown in Fig. 1141 The two distribu- 
tions agree well within the indicated errors. 



6.0.1 Comparison with standard chi statistics 

We also performed an exercise of applying the theory ly- 
ing behind the proper x^ distribution. We used the best 
fit value of Xdist{obs) from the Table Inland determined the 
error contour by adding the standard value of 4.61 (Avni 
1976) for 90% error, two parameters of interest (see Fig. llSH . 
Such a contour is much more compact than the contour de- 
termined directly from the simulations, as described above. 
This is surprising since the histogram of Xdist{i) approxi- 
mates quite well the theoretical x^ histogram (see Fig. 1161 
for the appropriate number of degrees of freedom (22 in the 
case of X-ray data). Error contours determined from simu- 
lations correspond to adding a value of order of ~ 8 — 28 
to the best fit Xdist{obs) . This shows that the Monte Carlo 
simulations are indeed extremely important, as argued by 
Uttley et al. (2002). 



APPENDIX B: STRUCTURE FUNCTION - 
ERROR DETERMINATION 

We again use two methods in order to determine the errors. 
The first method is to assign directly the error of the struc- 
ture function obtained observationally. The second method 
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Figure 15. Error contour at 90% confidence level for the bro- 
ken power law model of X-ray and optical NPSD determined as 
X^j„ -1-4.61; to be compared with thin lines in Fig.|^ upper panel. 
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Figure 16. The distribution of the Xdist(*) values (histogram) 
for 500 Monte Carlo curves in comparison with the theoretical 
plot for 22 dof. Parameters: X-ray data, model of a single break 
power law with the best fit parameters: log/o = —2.1, a = 1.5. 



includes the Monte Carlo simulations. Both methods are 
used in the present paper. 

B.l The direct method 

The structure function is computed in the timescale range 
determined by the data. We use logarithmic time bins with 
a width of 0.16. We determine the SF for bins with the 
number of pairs available greater than 7 pairs. The errors are 
assumed to come both from the flux measurement error and 



from the limited statistics of pairs. The flux measurement 
error is small, particularly for the photoelectric optical data 
(0.01, 0.015 and 0.02 - 0.025 mag for V, B and U band, 
correspondingly). The error on the data determined from 
the old photographic plates is below 20% and the X-ray 
data have errors of order of 10%. Therefore, statistical errors 
dominate. 

The statistical uncertainty in the SF(rfc) is equal to the 
ratio of the dispersion in the SF at a given bin to the number 
of pairs, i.e. 



' SF 



(st) 



K/2)2 



(19) 



where Xi 



is the mean pair difference and Uk is the num- 



ber of pairs in bin k. This error is plotted in Fig. |S| 

In our error analysis we included two additional effects, 
which should render the error determination more robust. 

For completeness, we took into account the errors asso- 
ciated with the flux measurements. This can be easy assessed 
through Monte Carlo simulations. Assuming that the errors 
in fluxes are normally distributed, we can modify each flux 
by random Gaussian deviates based on the quoted error for 
each data point. In a single Monte Carlo realization each 
data point is modifled, and the SF is computed. Computing 
a large number of independent realizations (500 - 1000) we 
obtain the mean SF and the rms deviation of the mean SF. 
We call this rms deviation to result from flux randomization 
a^p '•^''(rfe) (see Peterson, Wanders and Horne 1998). The 
third source of errors is associated with the observational 
sampling of the light curves. The specific times of observa- 
tions are in a sense 'randomly distributed' and the SF can 
be less sensitive to the choice of individual data points than 
the power spectrum. We can estimate these uncertainties 
through considering subsets of the original "parent" data 
set. This method is known as the bootstrap method. It can 
be used to evaluate the significance of the SF itself, its rms 
and other parameters such as slope " 6" of SF. The method 
works as follows: 

We have a set of observations (tl,xl), (t2,x2),..., 
(tn,xn). We can take from this set n pairs of randomly se- 
lected points without regard to whether or not they have 
been previously selected. Thus we select a subset of the orig- 
inal data points, in which some data pairs repeat more than 
once and some of the original pairs are ignored. Next the 
multiple pairs are removed so the size of the selected sample 
is reduced. For this data set we again calculate the SF and 
we repeat this process many times (500 - 1000). Multiple 
realizations lead to a mean and standard deviation for the 
structure function based on a randomly chosen subset of the 
original data points. We mark this rms value describing the 
effect of random subset selection as agp ^""{rk). 

Our computations show that bootstrapping gives the 
most conservative values of uncertainties in comparison with 
the other two described above. This is probably due to the 
fact that the considered subsets have smaller number of 
points than the real light curve. 

The total uncertainty of the SF can be obtained if we 
assume that the errors add in quadrature. 

CrSF{Tk) = (tIp '•"'''(rfc) +(7lp ^''{Tk) +0%p ''""{Tk) . (20) 

These error bars are used in Figs.|H]and |U] This proce- 
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dure allows us also to estimate the best-fit value and uncer- 
tainties of such parameter as a slope 6 of linear part of the 
SF. 

We note here that all those uncertainties are model in- 
dependent. 



B.2 Monte Carlo simulations of light curves 

In order to estimate the errors of the SF without a significant 
loss of information appearing in the bootstrap method we 
also perform the Monte Carlo simulations of the light curve 
itself based on the specific SF model. 

The SF with the linear part having the slope 6 can 
be obtained within the frame of the Poissonian model also 
called "shot noise" model. It means that in a physical sys- 
tem the variations are due to the stochastic superposition 
of independent flares of a given basic shape randomly dis- 
tributed in time, but of various durations and amplitudes. 
Under certain conditions, imposed on the flare, the ampli- 
tudes and durations of such process may result in the first- 
order structure function SF{t) being a power law of r with 
an index b which takes on values from to 2 (e.g. Terebizh 
et al. 1989; Sergeev 1999, Cid Fernandes, Sodr'e, & Vieira 
da Silva 2000). 

Sergeev (1999) showed that if the number n of flares of 
duration T/ is given by a power law distribution, i.e. 



(21) 



and if the flare amplitude A(Tf) also depends on the flare 
duration as a power law. 



A{Tf) ~ T: 



(22) 



then the dependence of the structure function on the interval 
T asymptotically converges to a power law shape 



SF{r) 



b = + 2(3 + 2. 



(23) 



In the considerations of Sergeev (1999), the extension of the 
power law dependencies of the flare parameters given by 
Eauations l21l and l22l were limited by a natural physical con- 
straints: the fiare duration was allowed to vary only between 
Tf and Tp""" . The adopted shape of the fiares was Gaus- 
sian, but it was shown that the profiles of fiares are of no 
importance. Flares can have any profile: square, triangle or 
Gaussian (e.g. Sergeev 1999, Cid Fernandes et al. 2000). 

We performed Monte Carlo simulations using the soft- 
ware of Sergeev with 5 input free parameters: 

(i) dT(fiare)= the mean time interval between two fiares 

(ii) a = index in the dependence of fiare numbers on du- 
ration, 

(iii) /3 = index in the dependence of the fiare amplitude 
on duration, 

(iv) Tmin ~ minimal duration of flares 

(v) Tmax = maximal duration of flares. 

Using Sergeev's program for simulation of light curves, 
we modeled such light curves using various values of a and 
/3. It is important to stress, however, that it is impossible 
to obtain separately the values of the parameters a and /3. 
Therefore, we took into account the relation b = 2 + a + 2(5, 
where 6 is a slope of the SF determined from observations. 

Changing all 5 parameters we simulated a number of 



light curves (N~500 - 1000) in each energy band and we 
computed the SF for each simulated light curve in the same 
manner as we did for the observed light curve. From each bin 
of all calculated SFs we subtracted 2Derr with D^rr being 
the mean variance of the measurement uncertainties. The 
model SF (corresponding to the given model and sampling 
pattern of the observed light curve) is then given by the 
mean of the N simulated SFs, SF{T)sim- The model error is 
given by the rms spread of the simulated model population 
around the mean value at each time interval r, (T(r)sim. 

The quality of the model representation of the observa- 
tional light curves is estimated with help of statistics Xtest) 
which we calculate as follows: 



Tf 



2 



SF{T)obs 



(24) 



We assumed that the best parameters of model are 
those yielding the minimum value of Xteat and the values 
of these parameters are listed in Table 0] 

This paper has been processed by the authors using the 
Blackwell Scientific Publications style file. 
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